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Abstract 

Two-dimensional dielectric microcavities are of widespread use in microoptics applica- 
tions. Recently, a trace formula has been established for dielectric cavities which relates 
their resonance spectrum to the periodic rays inside the cavity. In the present paper we 
extend this trace formula to a dielectric disk with a small scatterer. This system has been 
introduced for microlaser applications, because it has long-lived resonances with strongly 
directional far field. We show that its resonance spectrum contains signatures not only 
of periodic rays, but also of diffractive rays that occur in Keller's geometrical theory of 
diffraction. We compare our results with those for a closed cavity with Dirichlet boundary 
conditions. 
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1 Introduction 

Two-dimensional optical microcavities are flat dielectric objects that are surrounded by a 
medium of lower refractive index. They trap light by total internal reflection and are of great 
interest as miniature lasing devices with many technological applications |Vah04] . The con- 
finement of light results in long-lived states which is key for achieving low lasing threshold. In 
addition to the light confinement another highly desirable design feature is the directionality of 
the output. Examples of recent studies to achieve these goals include deformations of the bound- 
ary |DBD+08l lSHW+09j or insertion of a hole j WHOGj or a defect |AR04l ITVOTl [DM5W08] . A 
more detailed list of references can be found in these articles. 

The emission patterns of dielectric cavities are determined by the eigenmodes that are se- 
lected during the pumping process. Because dielectric cavities are open systems, the eigenmodes 
are resonance states with complex wavenumbers k = kr + iki. The real part of the wavenumber 
is related to the frequency u hj k^ = oo/c, where c is the speed of light, and the imaginary part 
is related to the lifetime r by /cj = — l/2cr. In practice, the resonance spectrum of a cavity has 
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Figure 1: Circular dielectric cavity of radius R and refractive index n in an outer medium 
of refractive index n. For the closed system we impose Diriclilet boundary conditions at the 
interface. For the perturbed systems a point scatterer is added on the x-axis at a distance d 
from the centre of the cavity. Throughout this article we will take n = 1. 

to be determined numerically or experimentally, because the only analytically solvable systems 
are the integrable cases with rotational symmetry like the circle. One interesting question is 
how one can understand the distribution of resonances in the complex plane that are found 
by the numerical or experimental methods. For this purpose, it has been very useful to adapt 
short-wavelength approximations that are of widespread use in the field of quantum or wave 
chaos |Sto99] . In a recent paper short-wave methods have been applied to determine a Weyl law 
for the number of resonances, a lower bound for their imaginary part, and to reveal signatures 
of periodic rays in the resonance spectrum |BDS08j . 

In this article we will focus on the trace formula for dielectric cavities that was obtained in 
IBDSnSj . It considers the excess density of states which is a sum over Lorentzian terms for the 
resonances. In the short-wave approximation this quantity is approximated by a smooth Weyl 
term, describing the mean density of resonances, plus a sum over the periodic rays inside the 
cavity. So far its main application has been to consider Fourier transforms of resonance spectra 
and identify the positions of its peaks with lengths of periodic orbits |BBD"'"10[ |BDD"'"09] . 



Trace formulas for dielectric cavities are on a less firm theoretical footing than the re- 
lated trace formulas for closed systems that were developed in the 1970's by Gutzwiller and 
Bahan/Bloch |BB7n[ iBBTll [BB721 [GutTll |Gut90] . There is not yet a systematic way of deriv- 
ing it, except for the integrable circular disk. It is known to require corrections for periodic rays 
with reflection angles that are close to the critical angle of total reflection. This leads to large 
finite size effects if only a small part of the resonance spectrum is known. A test of the true 
short-wave regime requires a large number of resonances which are difficult to obtain except 
for the circular disk. 

In this article we will study the trace formula for a dielectric circular disk with a point 
scatterer, see Fig. |ll It was shown in |DMSW08l IDMSWOQa] that this system has resonance 
modes with highly directional far field in addition to large quality factors Q = —kr/2ki. For our 
study it has the advantage that it is non-integrable but allows, nevertheless, the determination 
of a large number of resonances by a Green function method. 

The addition of a point scatterer leads to additional contributions in the trace formula from 
diffractive rays which are closed orbits that begin and end at the scatterer. These rays appear 
in Keller's geometrical theory of diffraction |Kel62] . They have been studied in detail for closed 
billiards |VWR,94l lPS95l ISiM IRRFOSi IRFM ILLMOBj . and we adopt these results here for 
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open dielectric systems. The closed orbits in a circular cavity have been recently classified in 
|BR09j ■ They undergo bifurcations as the position of the scatterer is changed which we will 
treat by uniform approximations. We will also consider the dielectric circle without scatterer, 
and we compare our results to those for closed cavities with Dirichlet boundary conditions 
(mirror boundary). 

Our article is organized as follows. After an introduction of the trace formula in section |2] 
we first consider the closed cavity without scatterer in section |3} and with scatterer in section |4j 
We then investigate the open dielectric cavity without scatterer in section |5| and with scatterer 
in section |6l Section [7] contains our conclusions. 



2 Trace formulas: from closed to open cavitities 

A common approach for studying the emission properties of microresonators is to consider 
them as passive systems and neglect the interaction between the electromagnetic field and the 
cavity material. The electromagnetic field for these dielectric cavities is then described by the 
three-dimensional vectorial Maxwell equations. 

Maxwell's equations simplify considerably if the electromagnetic field does not depend on 
one spacial direction, say in systems with cylindrical symmetry |Jac98j . The electromagnetic 
field can then be decomposed into two independent polarizations, called TM and TE modes. 
For each of these polarizations. Maxwell's equations can be reduced to a two-dimensional scalar 
equation in the xy-plane. For homogeneous dielectric cavities with refractive index n and cross- 
section D, this equation has the form, see e.g. |DMSW0"9a] . 

{V' + n\r) ij{r) = , n{r) = h ' ^ ^' (1) 

where r = (x, y) and k is the wavenumber. Here it is assumed that the surrounding medium 
has refractive index 1 (like air). The function ip{r) describes the z-component of the electric 
field in the case of TM modes, and the 2;-component of the magnetic field in the case of TE 
modes. All other components of the electromagnetic field can be obtained from the solutions 
of ([T]). The two modes satisfy different connection conditions at the boundary of the cavity. In 
this article we will only consider TM modes for which ip{r) and its derivatives are continuous 
across the boundary of the domain D. 

The reduction of Maxwell's equations to equation ([T]) requires cylindrical symmetry, but it 
is also approximately valid in thin dielectric cavities in which the height is of the order of the 
wavelength. In this case n is an effective refractive index that depends on the material as well 
as the thickness of the cavity. Although theoretical models exist for this index |LDA"'"07] . it 
is usually determined experimentally. Errors of this two-dimensional approximation have been 
investigated in |BDM0+n9j . 

Equation ([T]) is the basic equation that is widely used for the study of dielectric cavities. It 
describes open systems in which the wavefunctions extend to infinity. A related closed system 
is given by the Helmholtz equation for waves in a bounded domain D 

(V^ + A;^) ^(r) =0, r e D, (2) 

with Dirichlet boundary conditions for which the wavefunction vanishes at the boundary. This 
models, for example, an optical cavity with refractive index n = 1 and a mirror boundary. 
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Alternatively, it also agrees with the time-independent Schrodinger equation for quantum bil- 
liards. These are quantum versions of systems in which a point particle moves freely inside the 
domain D and is specularly reflected at the boundary. Quantum billiards have been studied 
extensively, and asymptotic methods for the semiclassical (short-wavelength) limit — )■ oo are 
well established. These methods have been very helpful in obtaining related short-wavelength 
approximations for dielectric cavities. For this reason we will discuss trace formulas for quantum 
billiards before we return to open dielectric cavities. 

The density of states is a useful characterisation of the eigenvalue spectrum of the Helmholtz 
equation and it is given by 

oo 

d{k) = ^ - , k> 0, (3) 

m=l 

where km are the positive fc-values for which equation ^ with Dirichlet boundary conditions 
has non-trivial solutions. The density of states is formally related to the Green function of the 

system by 

d{k) = -—lmtiG{r,r',k), (4) 

TT 

where the Green function is the solution of 

(V^ + k^) G{r,r',k) = 6{r-r'), r,r' e D, (5) 

with Dirichlet boundary conditions. In the limit where the wavelength is short compared to the 
dimensions of the billiard the Green function can be approximated by a sum over all classical 
trajectories connecting r' to r in the billiard system. If one then evaluates the trace in Q in the 
asymptotic regime of large wavenumbers one finds that the density of states is approximated by 
a sum of a smooth part do{k) and an oscillatory part dosdk) |BB70t lBB71t IBB72t IGut7T| IGutQOj 

d{k) ^ do{k) + dosc{k) . (6) 

The smooth part comes from the contributions of orbits of zero length. It is given in terms of 
the area A and perimeter L of the billiard system by Weyl's law |BH76j 

doik) = ^-^ + oil). (7) 

The oscillatory part consists of a sum over the classical periodic orbits in the billiard system. 
In this article we will consider also perturbations of quantum billiards and dielectric cavities 
by a point-scatterer inside the domain D. In this case there are additional contributions to 
the oscillatory part from so-called diffractive orbits. These are trajectories which start and end 
at the scatterer. Between start and end point they can have multiple reencounters with the 
scatterer, and at each encounter the outgoing and incoming angles are unrelated to each other. 
This can be related to the fact that an s-wave scatterer scatters uniformly in all directions. We 
hence write the oscillatory part as 

dosc{k) = dpo{k) + ddo{k) , (8) 

where dpo{k) gives the contribution of periodic orbits and ddo{k) gives the contribution of 
diffractive orbits. For the unperturbed cases (without scatterer) we have that ddo{k) = 0. Both 
periodic and diffractive trajectories give oscillatory terms to the density of states of the form 

d^ik)=x^A^ik)e''^^+c.c, (9) 
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where is the length of the trajectory ^, is an amphtude factor, and c.c. denotes the 
complex conjugate. Furthermore, Xs. is a phase factor that arises from the Dirichlet boundary 
conditions. It equals (—1)''^ where is the total number of reflections at the boundary. The 
contribution of a diffractive trajectory is weaker than that of a periodic one. The amplitude 
contains a factor for every scattering at the scatterer. 

In |BDS08j it was investigated how one can generalize the trace formula to dielectric cavities. 
This is not obvious, because dielectric cavities are open systems and equation ([T| admits a 
continuous spectrum corresponding to the scattering solutions. Using the Krein formula |Kre53t 
IKre62] it was shown that the relevant quantities to consider are the resonance states. These 
are solutions of equation ([T| with outgoing boundary conditions 

V^(r) oc e'^^l'l, |r| ^ cx). (10) 

They exist only for a discrete set of complex values of k with negative imaginary part. These 
complex fc-values coincide with the poles of the scattering matrix, and also with the poles of the 
Green function, and we denote them again by km- It was shown in |BDS08j that for dielectric 
cavities, one can again obtain a trace formula of the form ^ if one replaces the delta-functions 
on the left-hand side by Lorentzians 

^(^) = - - |Z (A; _ Rek^y + (Imkm)^ ■ ^^^^ 

Km. 

There is, however, an issue concerning which resonances to include in this sum. For TM modes 
one can typically distinguish two kind of resonances whose wavenumbers are well separated 
in the complex plane. The inner or Feshbach resonances have a high concentration inside 
the cavity and their wavenumbers lie in a strip —7 < Imkm < where 7 depends on the 
cavity. The outer or shape resonances are mainly concentrated outside the cavity and their 
wavenumbers have larger negative imaginary parts. Hence the outer resonances give only a 
smooth contribution to the density of states. These two types of resonances differ also in their 
behaviour as the refractive index n — )■ 00 as will be discussed later. 

In |BDS08j the smooth density of states do{k) on the right-hand side of the trace formula ^ 



was derived for the inner resonances. Hence one should also include in the sum (11) only the 
inner resonances. The result of |BDS08j for dQ{k) was inferred from the example of a circular 
dielectric cavity which is integrable and is given by 

do{k) = ^— + r{n)—, (12) 
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where the function r(n) has the form 



r{n) = l + - T^r^dt--/ -^dt (13) 



with 



The oscillatory part in the trace formula is determined by the ray dynamics inside the cavity. 
Since the contributing objects are formally the same as the trajectories in billiards we will adopt 
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the same terminology and speak of trajectories or orbits instead of rays. In this language, the 
contributing trajectories are exactly the same as in the closed cavity. The contribution of any 
single trajectory is modified from (|9| to 

d^{k) =nR^A^{nk)e'''''^^ + C.C.. (15) 

Apart from the additional factor given by the refractive index n, the main difference to the 
closed case is that the phase factor is replaced by the quantity which is the product of 
the Fresnel reflection coefficients for all the reflections of the trajectory at the boundary. If 
periodic orbits come in continuous families, then one should replace by an average (R^) over 
the family. For comparison with |BDS08j . note that we consider d{k) instead of d{E) = d{k)/2k. 

The trace formula for dielectric cavities is on a less firm basis than that for quantum billiards. 
The trajectory contributions ^ have been derived only for the case of the integrable circular 
cavity, and they have been generalized to other cavity shapes by physical considerations. It 
is also known that the contribution ^ becomes inaccurate if one of the reflection angles of a 
trajectory is close to the critical angle of total reflection |Jac98j 

9,{n) = sin-^ (^^^ , (16) 

where the reflection angle is defined with respect to the normal to the boundary. This inaccuracy 
is due to the Goos-Hanchen effect and Fresnel filtering |BDS08] . 

For TE modes, the separation of resonances into inner and outer resonances is less obvious 
than for TM modes |DBD"'"08[ IDMSW08] , and a trace formula has not yet been derived for TE 



modes. 

We note that for all our numerical computations in this paper we choose the radius R of 
the cavity to be 1. This is no restriction as the Equations ([T]) and ^ (together with their 
boundary conditions) scale with R. 



3 The closed circular cavity 

Short wavelength approximations are much better established for closed cavities than for open 
dielectric cavities. We consider in the next two sections a circular disk with Dirichlet boundary 
conditions, first without and then with perturbation. This allows us to introduce asymptotic 
short-wave approximations in a well-established setting and it will illustrate the accuracy of 
the approximations in these cases. It will set a benchmark for comparison with results for the 
dielectric disk discussed in Sections |5] and [6l 

The Helmholtz equation ^ for a circular disk with radius R is separable in polar coordi- 
nates, and we label the (unnormalized) solutions as 

i^m,q{r, 0) = Jm{km,qr) { ^^^^^^^j ' , ^ ^ ' (17) 

I cos(m0j , +m G i^o 

where Jm are Bessel functions of the first kindj^ The Dirichlet boundary conditions lead to the 
eigenvalues fc^.g = jm,q/R, where jm,q is the q^^ zero of Jm, q E N. Because the zeros of the 



^We choose this labeUng for convenience. The reader should not be confused by the fact that in the quantum 



mechanical interpretation the solutions ( 17) are not eigenfunctions of the angular momentum operator, but only 



of the squared angular momentum operator whose eigenvalues do not depend on the sign of m. 
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Bessel functions do not depend on the sign of m we have that the eigenvalues with m 7^ are 
doubly degenerate whereas those with m = are not degenerate. Hence the eigenvalues of the 
Helmholtz equation with Dirichlet boundary conditions are given by the multiset 

(T = {k^,q = jm,g/R : m G Z, g G N} . (18) 

In the next section we will also need the Green function for the circular disk. It can be 
obtained by taking the free Green function for the infinite plane and adding to it solutions 
of the Helmholtz equation ^ such that the Dirichlet boundary conditions are satisfied. The 
result is |SW71] 



Gir,r',k) = -J e'"^^^"'^V^(A;r 



JmikR) 



Hmikry) - /, Jmikr^^ 



= -'^Ho{k\r-r'\)+'- ^'"^''-'''jJ^MkrUUkr'), (19) 

m=—oo ^ ' 

where here and in the following Hm denotes the Hankel function -ffm"*, and r< and r> are the 
smaller and greater, respectively, of r and r' . 

As discussed in the previous section the trace formula ([6]) has the form 

d{k) = ^ 5(fc - fc^,,) ^ d^{k) + rfosc(fc) , (20) 

where the sum is taken over the eigenvalues as described above. We will apply the trace formula 
to investigate signatures of periodic orbits in the spectrum of the cavity. For this purpose it 
is useful to consider the length spectrum which is given by the Fourier transform F{1) of the 
density of states. Since the trajectories give oscillatory contributions of the form ^ to the 
trace formula, F{1) shows a peak structure at periodic orbits. In practice we can numerically 
determine the spectrum only within a finite window < A; < fcmax- With this available range 
in mind, the transform is weighted with a Gaussian cut-off function W{k) in order to reduce 
the influence of larger eigenvalues. We choose W{k) = exp(— tfc^) with cut-off parameter 
t = 10 /k^^^. In our numerical evaluations fcmax = 100 which corresponds to a window with 
2456 eigenvalues. 

Since d{k) is a comb of delta functions its weighted Fourier transform is 

poo 

F{1) = / d{k) W{k) e-^'^' dk = y] g-''^^.^'-'^-..* , (21) 
^0 

in,q 

We will compare this result to its short-wave approximation which is obtained by evaluating 
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the weighted Fourier transform of the right-hand side of the trace formula ( 20 ) . 

The smooth part of the density of states do{k) is given by equation ([7| with area A = ttR 
and perimeter L = 2'kR. For brevity we will only consider the real part of the length spectrum. 
The real part of the weighted Fourier transform of the smooth part is then given by 



where erf(-) is the error function. 



A e-''/4* lAl J il 
L erf 



t \2Vi 



(22) 
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The oscillatory part of the density of states can be obtained by performing the trace in ^ 
for the Green function of the circular disk (19) asymptotically as A; — )■ oo. It is given by a sum 



over periodic orbits and has the form (see e.g. |PVC01j ) 



2 °' ^'^'^ [2k 
dpo{k) = -J2J23pA\ -r cos {kip - Up7i/2 + 7r/4) (-l)'' , (23) 

^ r=2 w=l V ^ P 

where [-J is the floor function (integer part). The periodic orbits in the circular disk occur 
in families because of the rotational symmetry and we use the integers r and w to label the 
different families. For simplicity we combine them into a single index p = {r,w). The orbits 
have w rotations around the billiard centre and r reflections at the boundary. For w = 1 these 
are regular polygons with r vertices. In Fig. [2] we show representative periodic orbits for some 
of the families with w = 1. Each periodic orbit PO{r,w) reflects from the boundary with an 
angle of incidence 6p measured from the normal to the boundary, given by 

We occasionally denote it by 9r,w The periodic orbit length is 

lp = 2rRcos{ep) , (25) 
and the area of the annulus filled by all the members of each family in configuration space is 

Ap = nR'^cos'^{ep). (26) 



The Maslov index Up in (23) counts the number of conjugate points and is simply given by 
Pp = r. The factor of (—1)^ comes from the Dirichlet reflections at the boundary as discussed 
earlier. Finally Qp is a degeneracy factor due to time-reversed orbits. We have Qp = 1 for the 
diameter orbit and its repetitions and Qp = 2 for all other orbits. We note that the sum over r 



and w in (23 ) can be decomposed into a sum over relatively prime r and w and integer multiples 
of the two. The integer multiples then correspond to repetitions of a so called primitive periodic 
orbit characterized by the relatively prime r and w. 

The contribution of the periodic orbits to the length spectrum is obtained by taking the 



Fourier transform with the weight function W{k) of (23). The integral can be performed 
analytically and its real part is given by 

1 A (-lU(3'-+i)/2J /2\3/4 



P 

X ~' 



{Ip + l), 



where denotes the parabolic cylinder function. In Fig. |3} we plot the oscillating part of the 
semiclassical length spectrum, ReFosc(/), together with the exact length spectrum minus the 
smooth part Re[F(/) — Fo{l)]. One sees that the short-wave approximation describes the exact 
curve very well, and it shows periodic orbits as peak structures at their corresponding lengths. 
Note that for fixed rotation number w we have Ip — )■ 27rw as r — )■ oo, and hence the periodic 
orbits approach the limit of uniform circular motion around the inner edge of the billiard. The 
orbits therefore accumulate as one approaches the values / = 2w7r from below in the length 
spectrum, and therefore the orbits PO(r, 1) cannot all be resolved in Fig. [3] 
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Figure 2: The four shortest periodic orbits in the circular cavity: the diameter, the equilateral 
triangle, the square and the regular pentagon. All of the orbits occur in continuous families 
obtained by rotation about the origin. The angle 9^^^ is shown for the triangular periodic orbit. 




Figure 3: Length spectra for the unperturbed circular billiard. The exact spectrum is computed 
from (21) with unperturbed eigenvalues (minus the smooth part (22)) and the short-wave 
approximation from (27) with periodic orbits {t = 0.001). 
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4 The closed cavity with point scatterer 



4.1 Exact solutions 

We perturb the system in the previous section by inserting a point scatterer. Without loss of 
generahty we locate it on the x-axis at a position d = (d, 0) in Cartesian coordinates, see Fig. [l] 
Self-adjoint extension theory |Zor80] tells us that the perturbed eigenvalues of this system are 
solutions of |Shi94 

1 - AG,cg(d,d,fc) = 0, (28) 

where A is a coupling strength parameter and Gregld, d, fc) is the regularised Green function. 
The regularized Green function can be obtained from G in ((T9| by subtracting its divergent 



part. Replacing the Hankel function Hq in (19) by its asymptotic representation as r' — ?■ r = d 



and subtracting the logarithmically divergent term gives 

G„.(d, d, = -i + i (log 4 + 7) + ^ E , (29) 

^ m=—oo ^ ' 

where 7 is Euler's constant and fco is an arbitrary constant corresponding to a choice of the 



regularisation. The imaginary part of (29) vanishes, and we find that the perturbed eigenvalues, 



/c, are determined as the solutions of /(/c) = 0, where / is the transcendental function 

(,og^+.)4|;^.„.jM (30) 



with 




(31) 



The parameter a in (30) is a coupling constant which is a measure of the strength of 
the perturbation. It replaces the previous parameters A and /cq and is related to them by 
A = — 27r/ In(fcoo)- It has a direct physical interpretation in the limit a — )■ where the system 
with the point scatterer approaches that with a small scattering disk of radius a with Dirichlet 
boundary conditions jRWW96| . Further ways to model the point scatterer experimentally will 
be discussed in section 16.11 The value of a ranges from to 00 which correspond to A — )■ 0+ 



and A — 7- 0„, respectively. At these two limits the zeros of ( |30[ ) approach the unperturbed 
eigenvalues of Section |3| 



The perturbation described by (28) is a rank one perturbation. In the case of degenerate 
eigenvalues of the original problem this means that only one of the eigenvalues of a degenerate 
set changes due to the perturbation while all others remain unaffected. For the unperturbed 
circle all eigenvalues with (m 7^ 0) are two-fold degenerate, and only one of each pair changes 



according to (30) whereas the other remains as it is. This can also be seen from the solutions 
(17), because the wavefunctions with angular dependence sin(m</)) have a nodal line on the 
X-axis and are unaffected by the presence of the scatterer. The spectrum a of all eigenvalues 
of the perturbed system is therefore 

o = {km,q = km,g G (T : -TTi, g G N} U {krr,^q G M : m G Nq, ? G N, fCkm,q) = 0} , (32) 

where we label the km,q as follows. We give them the same m and q numbers as the unperturbed 
level at which they arrive as a —t- 00. The reason for choosing the a — > 00 and not the a — )■ 
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Figure 4: Level dynamics of the perturbed eigenvalues of the circular cavity with point scatterer 
positioned at ci = 0.59 for k G [0, 12]. The solid curves are the eigenvalues which change and 
the dashed horizontal lines are the eigenvalues which do not change as a is varied. 



limit for the labelling is that otherwise we would miss the first perturbed level. This eigenvalue 
starts at = fco,i at a = oo, and it decreases with decreasing a and finally becomes complex. 
A similar effect occurs also in one dimension where a delta function potential with negative 
coefficient supports exactly one bound state |Lev02j . 

As an illustration we plot in Fig. |4] the eigenvalues of a circular disk with point scatterer 
at d = 0.59 as a function of the parameter a. They have the property that there is always 
one perturbed eigenvalue between any two adjacent unperturbed ones. This follows from the 
fact that the Green function is a monotonically decreasing function between its poles. We used 
this property to solve (28) using a bisection method. Due to the logarithmic dependence on 
the strength in (30), we vary a exponentially from a = 10^^ to a = 10"^". Note that the first 
perturbed eigenvalue is complex beyond a certain value of a. 

We consider now the length spectrum ( [2l] ) of the perturbed cavity. It is of the same form 
as before, namely 

= ^e-^'^™-'-*^^-% (33) 

m,q 

where the sum is now over the eigenvalues of the perturbed disk. The parameter t is again 
t = 10/ kl^^^ with fcmax = 100. In Fig. p\ we plot the unperturbed and perturbed length 



spectra together, with the smooth part Fo{l) in (22) subtracted. We see that the periodic 
orbit structures are the same for both curves, but there also are small discrepancies. We will 
see in the next section that these discrepancies can be attributed to diffractive orbits. In 
order to visualize their contribution better we will consider in the following the difference, 
AF(/) = F{1) — F{1), of the perturbed and unperturbed length spectra. 
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Figure 5: Exact length spectra for the unperturbed and perturbed circular cavities. The 
unperturbed spectrum is computed from (21) and the perturbed spectrum from (33). Both 
spectra have had the smooth part subtracted {d = 0.59, = 1.0, t = 0.001). 



4.2 Short-wave approximation 

The short wavelength approximation for the perturbed system can be obtained again from its 
Green function. This Green function follows from self-adjoint extension theory as 



^(r, r', k) = G(r, r', k) + G(r, d, k) 



X 



l-AGrcg(d,d,A;) 



G(d,r',A;), 



(34) 



where G is the Green function of the unperturbed disk. If one then evaluates the trace of the 
Green function in ^ to obtain the density of states asymptotically as — )■ oo, one finds that 
the first term on the right hand side of (34) gives the same contributions as in the unperturbed 
case (mean part and periodic orbit terms). The second term gives additional contributions. 
If one expands the fraction in a geometric series one finds that these contributions come from 
trajectories that start and end at the scatterer and do this an arbitrary number of times. These 
are the diffractive trajectories that occur in Keller's geometrical theory of diffraction |Kel62] . 

For a single diffractive orbit, labelled by d, the contribution to the density of states is 
|VWR94j 



dd{k) = Re 



I. 



n 



Vik) 



\ v/8vrA;|(Mrf,,)i2| 



exp 



i kL„ 



_ 3/irf7r 
2 4 



V'd 



(35) 



where j labels the different parts of the diffractive orbit between consecutive encounters with 
the scatterer. The total length of the diffractive orbit is L^, the Maslov index counts the 
number of conjugate points, is the total number of reflections with the boundary and (Mj)i2 
is the 12-element of the monodromy matrix Mj. The total number of encounters of a diffractive 
orbit with the scatterer is fid- Furthermore, the length 1^ is the total length L^, divided by 
the repetition number in the case that the diffractive orbit is a multiple repetition of a shorter 
one. The diffraction coefficient is V, which for a point scatterer is independent of incoming and 
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a 



Figure 6: Point scatterer position c? as a function of the initial angles a of diffractive orbits 
DO{r,w). For fixed r and w, the function (39) is plotted (solid lines) for the diffractive orbits 
and the function (46) is plotted (dashed lines) for the periodic orbits. The radial diffractive 
orbits R0(2/c + l,k + 1) and R0(2fc + 1, k) correspond to the vertical lines a = and a = tt, 
respectively. 



outgoing angles and given by |RWW96| IEv96] 

V{k) 



27r 



i7r/2 - 7 - log(A;a/2) ' 



(36) 



Finally we should mention that if a trajectory goes straight through the scatterer without 
changing its angle, then it counts twice, once counting the point of coincidence with the scatterer 
as encounter and once not counting it as such. In our computations we chose a value of a = 1 
for our perturbation parameter. This value is suitable for seeing clear diffractive effects for 
wavenumbers in the range [0, 100]. 

We investigate now the different kinds of diffractive orbits that occur in the circular cavity. 
It is sufficient to consider diffractive orbits that have just one encounter with the scatterer, 
because all other diffractive orbits can be obtained from combinations of them. These single 
encounter diffractive orbits occur in two forms: they do or do not coincide with a periodic orbit 
that goes through the scatterer. We investigate first the latter type and generically label them 
as DO(r, w), where r is again the number of reflections at the boundary and w is the number of 
rotations around the billiard centre. Each diffractive orbit has r + 1 legs. This gives the angle 
of incidence 9 of diffractive orbits as solutions of the transcendental equation |BR09j 



/rvr 



r0 



(37) 



The simple dependence on the rotation number w in the condition (37) means that the diffrac- 



tive orbits are separated into two classes; one class consists of orbits with an odd number of 
rotations, the other of orbits with an even number. It is more convenient to write the condition 
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RO(1,0) 



R0(1,1) 



R0(2,l) 




Figure 7: The first few diffractive orbits by lengtli of tlie circular billiard with point scatterer. 
The scatter is positioned at ci = 0.59 in all panels except for D0(4, 1), where it is positioned 
at d ~0.6825, the bifurcation point where the real pair of orbits D0(4, 1) and D0(4, 1)' come 
into existence. For D0(3, 1), we also show the corresponding diffractive periodic orbit P0(3, 1), 
with its reflected partner (dashed lines), which comes into existence ai d = 0.5. 



in terms of the initial angle a at which each diffractive orbit starts with respect to the symmetry 
axis. The relation between this angle and the angle of incidence is 



9 = 



[r - Dtt 



iw — l)7r — a 



This gives the position of the scatterer as a function of the initial angle as 



d 



R 



sm 



sm a 



(r — l)7r 
2r 



a 
r 



(38) 



(39) 



In Fig. 6 we plot the condition (39 ) for a few orbits that we will later use. For a particular choice 



of r, w and d^ this condition is solved numerically for a by a Newton procedure. Depending on 
the choice of d^ we see that in some instances there is no solution for the diffractive orbit and 
in others there are two solutions. In this way, some orbits are created by tangent bifurcations 
as the parameter d is varied and we will discuss this later. 

There are two radial diffractive orbits in the circle billiard both travelling along the x-axis 
and both having r = 1 reflections with the boundary. The longer of these two orbits passes 
through the centre of the billiard and we hence label it as RO(r, w) = R0(1, 1). The shorter 
orbit does not pass through the centre and we hence label it as RO(r, vS) = R0(1, 0). Either of 
these two orbits can be added to multiples of the radial periodic orbit P0(2, 1) where we regard 
P0(2, 1) as always missing the scatterer. In this way, these combination orbits are primitive 
as they always only contain one encounter with the scatterer. We therefore denote all these 
orbits by R0(2/c + l,k + 1) and R0(2/c + 1, k) for k G No, for the longer and shorter orbits 
respectively. 

The first and last legs of each diffractive orbit are of equal length 

X = d cos a + R cos 9 . (40) 
The remaining r — 1 legs are of length / = 2Rcos9, giving the total length of the orbit as 

ld = 2{rRcos9 + dcosa). (41) 
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Using the basic monodromy matrices for motion on a straight path and reflection at a curved 
boundary [SS90] we find the 12-element of the monodromy matrix to be 

Mi2 = „^ ^ {Rcos9-x){{r-l)Rcos9-rx). (42) 
li cos 



The Maslov index in (35) is the number of conjugate points and can be computed from the 



number of times M12 changes sign along an orbit path. The result is 

r-2, x<{l-l/r)/2 
Vd={ r - 1 , (/ - //r)/2 < X < //2 , (43) 
r , x> 1/2 

which also remains valid for the radial orbits RO(1,0) and R0(1, 1) even though the length 
/ = 2R does not form part of these two orbits. Due to the time-reversal symmetry in the billiard 
system each diffractive orbit comes in a pair, except for the radial orbits RO of which there is 
only one of each. 

As mentioned before, diffractive orbits may also occur in the form of periodic orbits, and 
we denoted them also by PO(r, w), depending on whether the condition 

holds, where r is again the number of encounters with the boundary and w the number of 
rotations around the disk's centre. For a given periodic orbit family, for which each member 



has an angle of incidence at the boundary given by (24), the value of d at which we have equality 



in (44) is the caustic radius of the family. There are two isolated members of each family that 
serve as diffractive orbits which are related to each other by symmetry with respect to the x- 
axis. For both these orbits, the first and last legs have lengths Xi^2 = i ^/^^ + 4(6?^ — -R^))/2, 
respectively, where we again have / = 2Rcos6 and the 12-element of the monodromy matrix 
for diffractive-periodic orbits is given by 

= ^^JR' sin' e - d') . (45) 

The corresponding number of conjugate points is = r — 1. Each of the periodic orbits is 
doubly degenerate due to the reflection symmetry and in turn each of these is doubly generate 
due to time reversal symmetry. Each primitive periodic orbit must therefore be counted four 
times. 

Similarly to the diffractive orbits, we can give an expression for the position of the scatterer 
in terms of the initial angle a at which each periodic orbit starts. It is given by 

d = R'^. (46) 
sma 

We plot this condition in Fig. [6] alongside the corresponding curves for the diffractive orbits for 
a couple of orbits that we will use later. We see that the curves corresponding to the periodic 
orbits intersect certain curves corresponding to the diffractive orbits and in this way some orbits 
are created by pitchfork bifurcations as the parameter d is varied. We will give a review of all 
the bifurcations later. 
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We have so far only considered diffractive orbits which have one encounter with the scatterer. 
In this case we have /i^ = 1 and Id = in (35). Diffractive orbits consisting of more than 



one encounter with the scatterer are simply combinations of the single encounter orbits that we 
have discussed. It is important to distinguish between multiply diffractive radial orbits that are 
combinations of the primitive orbits R0(1, 0), R0(1, 1) and P0(2, 1) and the singly-diffractive 
orbits R0(2/c + 1, k) and R0(2fc + 1, 1 + k). The diffractive contribution to the trace formula 
for the circular billiard with point scatterer is finally 



^^do(A;) = 5^Ad(-ir'^exp|i(^A;L, 



^--^ ^1 ^+c.c. 



(47) 



where the amplitudes are given by 




9d,p{k) 



L v/87rA;|(M,,,; 



(48) 



121 



and the sum runs over all diffractive orbits described above, labelled by d. The product runs 
over all primitive diffractive orbits RO(r, it?), 'DO{r,w) and PO(r, it?). In the amplitudes (48) 



we have also included a degeneracy factor g^^j which takes into account the double degeneracy 
of non-radial diffractive orbits and the quadruple degeneracy of the non-radial periodic orbits, 
as discussed above, i.e. 




if (d, j) = R0(2A; + 1, k), R0(2fc + 1, 1 + A;), A; G No 
if (d,j) = DO(r,ti;), P0(2, 1) 
if (ci,j) = PO(r,ti;), r > 3 



(49) 



The presence of the scatterer also adds a small additional term to the smooth part (io(/c) 



of the density of states. It can be obtained from the Green function (34). The smooth part is 
then given by 

" .2 / r_ _ \ 21 



do{k) 



Ak 
27 



L 

Atx 



1 
2k 



TC 

T 



+ 7 + log 



ka 



(50) 



For convenience, we slightly change the definition of the length spectrum for the perturbed 
cavity and consider 

poo 

(51) 



F{1) = Re / A:^— /2rf(A:) W{k)e-'''^dk , 
Jo 

where we introduce /imax as the chosen maximum number of encounters of all the diffractive 
orbits with the scatterer to be included in the short wavelength approximation. This removes 
the singularity in the integrand when d{k) is replaced by its approximation. The weight function 
W{k) is the same as in Section [3] In order to see the effect of the diffractive orbits in the 
perturbed spectrum we then consider the difference between the perturbed and unperturbed 
length spectra AF{1) = F[l) — F{1). In the short wavelength approximation this quantity is 
given only by the diffractive orbit contribution and the extra contribution to the smooth part, 
because the periodic orbit terms and the area and perimeter terms in the smooth part cancel. 
We hence have 



AF(/) 



Jo 



2A; 



/ , ka\ 



n -1 



W{k)e-'''^dk. (52) 
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This integral will be computed numerically. In the lower panel of Fig. [5] we plot AF{1) for a 
scatterer position d = 0.59, scatterer strength a = 1 and /imax = 4 (full line) and compare it to 



the approximation (52) (dashed line). The diffractive orbits included are those given in Fig. [7| 
except D0(4, 1) and D0(4, 1)' which are not real at = 0.59 (we detail this below), plus the 
radial periodic orbit P0(2,l). We see that the peaks can be identified with the contributions 
of the diffractive orbits. The agreement is very good except at a few particular lengths. 

We have already mentioned that the diffractive orbits and isolated periodic orbits that 
serve as diffractive orbits are involved in pitchfork and tangent bifurcations. As with standard 



periodic orbit contributions, the diffractive contribution (47) to the trace breaks down near 
bifurcations of the diffractive orbits, because the matrix element M12 vanishes there. If one 
chooses a value of d in such a way that the scatterer is away from bifurcations then the length 



spectra difference AF is well described by the approximation (52). In order to deal with all 
values of d, one has to apply uniform approximations to the diffractive contribution to the 
density of states. The bifurcation scenarios for closed orbits (diffractive orbits) in a circular 
billiard have recently been described in |BR09] and we here summarise the results in Tab. [l} 

The bifurcations occur in three groups. Group I consists of diffractive orbits D0(2/c, k) (even 
r and maximum w) that are created from the break-up of the tori (BT) at d = oi the radial 
diffractive periodic orbits P0(2/c, k). Group II consists of diffractive orbits D0(2fc + 1, k) (odd 
r and maximum w) that are created by pitchfork (PF) bifurcation at d = R/r from the radial 
diffractive orbits R0(2/c+l, k). These orbits undergo a further PF bifurcation at d = Rcos{6r,w) 
to create the periodic orbits P0(2A; + 1, k) (note that this value of d is the caustic radius of 
the torus which contains the created periodic orbit). Group III consists of all other diffractive 
orbits. They are created in pairs DO(r, w) and DO(r, w)' by tangent (TG) bifurcation at dr^^, 
which corresponds to the value of each minimum in Fig. [6) One member of each of these pairs, 
DO(r, w), undergoes a PF bifurcation, again at the caustic radius d = i?cos(6'r,to), to create the 
periodic orbit PO(r, w). For group III, both orbits are not real before the bifurcation (leading 
to so called ghost orbits, see below). At the boundary, all the DOs become POs as indicated in 
the last column. Finally, we note that for the circle billiard there are no bifurcations involving 
only periodic orbits. 





bifurcation sequence group 


J r even 
w = r/2 


r odd 
• u; = (r-l)/2 


r even/odd 
u; < (r- l)/2 


1st bifurcation 


type 


BT 


PF 


TG 


at d = 





R/r 


dr,w 


from 


PO{r,w) 


R0{2k + l,k) 




to 


BO{r,w) 


BO{r,w) 


BO{r,wy 


BO{r,w) 


2nd bifurcation 


type 




PF 




PF 


at d = 




RcOs{6r,w) 




RcOs{9r^w) 


to 




PO{r,w) 




PO{r,w) 




at d = R 


PO{r + l,w) 


PO{r + l,w) 


PO(r - l,w) 


PO(r + l,tt;) 



Table 1: Summary of the diffractive orbit (DO) bifurcations for the circle billiard, which are 
separated into three groups. The table should be read columnwise and is explained in the text. 

Using the results in |SS97j for the uniform approximations of periodic orbits involved in 
pitchfork and tangent bifurcations, we can write down the corresponding formulae for the 
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Figure 8: Length spectra difference AF(/) for the perturbed circular bilhard. The exact 
difference (full lines) is computed from unperturbed and perturbed eigenvalues, and the short- 
wave approximation (dashed lines) from (52) with diffractive orbits. The upper panel in- 
cludes uniform approximations for the D0(3, 1) + P0(3, 1) pitchfork bifurcation (PF) and the 
D0(4, 1) -|- D0(4, 1)' tangent bifurcation (TG). These are labelled on the plot next to the cor- 
responding peak structure as are the other contributing orbits. The lower panel shows the 
short-wave approximation without uniform approximations, {d = 0.59, a = 1.0, /imax = 4, 
t = 0.001.) 
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uniform approximation of pairs of diffractive orbits involved in the bifurcations in Tab. [T| For 
our pitchfork bifurcations involving a diffractive orbit labelled by (exists before and after the 
bifurcation) and a pair of periodic orbits labelled by 1 (exist only after the bifurcation), the 
uniform approximation is 



dpp{k) = Re|27rA;AL|^/^(-l)"'*exp i (^kL - "^u - ^fi^ 



+ 



V2) 



(53) 



where L = (Li + Lo)/2 and AL = (Li — Lq)/2. The amplitudes Aq and Ai are given by (48), 
(Ti = sign(AL) and (72 = —1 before the bifurcation and (72 = 1 after the bifurcation. The 
index v is the average Maslov index of the orbits after the bifurcation. In our case we always 
have pitchfork bifurcations involving a diffractive orbit and a corresponding diffractive periodic 
orbit. 

For the uniform approximation of a tangent bifurcations, it is useful to complexify the phase 
space and view the birth of two (real) periodic orbits as to originate from two complex periodic 
orbits (so callled ghost orbits) prior to the bifurcation. Labeling the two real orbits by 1 and 2 
gives the contribution from the uniform approximation after the bifurcation 



d'VG{k) 



Re 



87rA;AL 



3 

A1 + A2 



— isign(AL) 



1/2 



1) exp 



(j_i/3 (A;|AL|) + Ji/3(A:|AL|)) 



A1-A2 



{j^2/3{k\AL\)-J2/s{k\AL\)) 



(54) 



where L and u are the average length and Maslov index, respectively, and AL = Li — L2. 



The ghost orbits before the bifurcation can be computed from the complex solutions of (39). 



It is then convenient to write the contribution from the uniform approximation before the 
bifurcation in the form 



dTG{k) 



Re 



8A;AL 



X 



TT 

A1 + A2 



1/2 



exp 



i ( /cL - % - ^^id 



K,„{k\AL\) + ^^^-^K2„{k\AL\) 



(55) 
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with L, AL and v as above and K being the modified Bessel function of the second kind. In 
this case, the lengths, actions and monodromy matrices are all complex valued. For the correct 
way to complexify, see |SS97] . 

In some cases we need to consider also multiply diffractive orbits in which one singly diffrac- 
tive part bifurcates. In these cases, the uniform approximations above still hold when inserting 
the quantities of the full multiply diffractive orbits. 

As the scatterer is moved from the centre of the billiard toward the edge, the bifurcations 
happen in sequence and, for a given value of one must determine which bifurcations are 
within its vicinity. In most cases this consists of determining which tangent bifurcation was 
last to occur and which is next to occur. For a fixed value of d we determine the next tangent 
bifurcation to occur in the direction of increasing d and input the corresponding ghost orbits 
into the diffractive orbit sum. 

In the upper panel of Fig. [8] we show the length spectrum difference AF together with the 
short-wave approximation using the uniform approximations discussed above for d = 0.59, a = 1 
and /imax = 4. Our choice of scatterer position is generic as it is not exactly at a bifurcation point 
and is between pitchfork and tangent bifurcations allowing us to apply both types of uniform 
approximation. Comparing this to the lower panel we see that the uniform approximations have 
improved the correspondence. In fact, we can now identify the discrepancies in the lower panel 
with particular orbits. We see discrepancies near the lengths corresponding to the diffractive 
orbit D0(3,l) (of length / = 5.2204), the diffractive periodic orbit P0(3, 1) (of length / = 
5.1962) that was created from it in a pitchfork bifurcation at d = 0.5, and the real part of the 
lengths of the ghost orbits D0(4, 1) and D0(4, 1)' (/ = 5.7196). For d = 0.5 the scatterer would 
be positioned exactly at the bifurcation point of the pitchfork bifurcation. In this instance one 
has to use the limiting form of the corresponding uniform approximation |SS97j . 

In Fig. [9] we give a three-dimensional plot of the exact length spectrum difference AF(/) for 



the scatterer varying from the centre of the billiard to the boundary, i.e. d G [0, 1]. In Fig. 10 
we plot the maxima of this plot in the {d, /)-plane together with the lengths of the diffractive 
orbits. We are able to identify all of the expected diffractive orbits. In particular, we can 
identify multiple diffractions of the shortest radial diffractive orbit RO(1,0) and combinations 
of longer orbits with this radial orbit. We can also identify the results of broken tori at d = 
and sequences of two pitchfork bifurcations whence the diffractive orbits D0{2k + l,k) with 
fc e No are created in the first bifurcation and the associated periodic orbits P0(2A; + l,k) 
are created in the second. Further, we can identify the sequences of tangent and pitchfork 
bifurcations whence the pair of diffractive orbits DO{r,w) and DO(r, w)' are created in the 
tangent bifurcation and the associated periodic orbit PO{r,w) is then created from DO(r, w)' 
in the pitchfork bifurcation. Additionally, for the tangent bifurcations, we see a trail of maxima 
leading from the bifurcation point that we identify with the ghost orbits included in the uniform 
approximation (55). 

In the TG-PF bifurcation sequences, the two bifurcations always lie very near to each other 
with respect to the bifurcation parameter d, which causes a problem when implementing the 
two corresponding uniform approximations. This has not yet been addressed and is remarked 
upon in |BR09] . As explained above, for a generic scatterer position we always consider the 
nearest bifurcation(s). For the value d = 0.59 used in the figures, this is pitchfork bifurcation 
at d = 0.5 and the tangent bifurcation at d = 0.6825. 
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Figure 9: Three dimensional plot of the absolute value of the exact length spectrum difference 
|AF(/)| for the perturbed circular billiard for varying position d of the scatterer, as computed 
from the resonance spectra, (a = 1.0, t = 0.001.) 
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Figure 10: Bifurcations in the circle billiard with point scatterer. The circles correspond to 
the maxima of the exact length spectrum difference in Fig. [9] and the lines are the lengths of 
the classical diffractive orbits. They are labelled as described in the text. Two bifurcation 
sequences are highlighted. The first sequence involves two pitchfork bifurcation: first the orbit 
D0(3, 1) splits from R0(3, 1) at PFi, and later P0(3, 1) splits from D0(3, 1) at PF2. In the 
second sequence DO (4, 1) and DO (4, 1)' are created in a tangent bifurcation at TG, and then 
P0(4, 1) splits from D0(4, 1) in the pitchfork bifurcation PF3. Before the tangent bifurcation 
we circle maxima of the length spectrum which correspond to ghost orbits. Note that for 
illustrative purposes not all classical orbits are marked in the plot. 
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5 The open dielectric cavity 



Having studied the closed cavity we now turn to the optical microcavity, which is a circular 
disk of dielectric material of refractive index n > 1 in an outside medium of lower refractive 
index n = 1. Waves can now reflect and refract at the interface and, restricting ourselves to TM 
polarisation, we consider the reduced Maxwell equations in the form of the Helmholtz equation 

0- 

5.1 Exact solutions 

The spectrum for outgoing boundary conditions consists of complex-valued resonances of the 
form k = kr + iki, with ki < 0. It can be found from the Green function. Inside the cavity, the 
outgoing Green function is given in polar coordinates r = (r, 0) by [DMSWOSt IDMSWOOa] 

i i °° C 

G{v,v' , k) = —-Ho{kn\r — r'\) + - cos [m((f) — (p')] Jrniknr)Jmiknr') . (56) 



m=— oo 



The coefficients in this expansion are given by 

Am = Jm{knK)H'^{kK) — nJ'^{knR) HmikR) , 

Cm = Hm{knR)H'^{kR) - nH'^{knR)Hm{kR) , (57) 

which, by using the formula for Bessel function derivatives jAS64j . we write as 

Am = nHm{kR)Jm+i{knR) - Jm{knR)Hm+i{kR) , 

Cm = nHm{kR)Hm+i{knR) - Hm{knR)Hm+i{kR) , (58) 



which is advantageous for numerical purposes. The spectrum is obtained from the poles of (56) 
which are given by Am = and hence each resonance is determined as a complex root of the 
function [DMSWOS] 

fmik) = nHm{kR)Jm+i{knR) - Jm{knR)Hm+i{kR) . (59) 

The corresponding eigenfunctions for a solution of fm{k) = 0, which represent the z-component 
of the electric field, have the form 

Hm{kR)Jm{knr) fsin(m0) , -m G N 

Jm{knR) lcos(m0), m G No ' ^ ^ 



for r < i? and 




^(r,0) = ff„(A:r)<;""^V^";' ^ , (61) 



for r > R. 

We mentioned in Section |2] that there are two kinds of resonances, inner and outer. One way 
to distinguish between them is to consider the limit as the refractive index n — )■ oo. In this limit 
the system becomes closed and separates into an inside and an outside problem. The outer 
resonances are those that correspond in this limit to solutions of the outside scattering problem 
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with Dirichlet boundary conditions at the disk |DBD+08l IBDS081 IDMSWOQbj . We will not 
consider them further, because they lie deeper in the complex plane than the inner resonances, 
and they do not play a role in the trace formula. The inner resonances on the other hand 
have the property that their wavenumbers become real in the limit n — oo. We label these 
resonances by two indices m and q where for fixed m, the index q enumerates the resonances 
in the order of increasing real part of their wavenumbers starting from q = 1. By considering 



the asymptotic behaviour of the function fm{k) in (59) one finds that the wavenumbers of the 



resonances labeled this way have the limiting behaviors |R,RP+n8[ IDMSWOQb] 



lim nkm,qR = jm-l,q, m^O 
n—>-co 

lim nko^gR = q^l (62) 

n— s>oo 

lim nko^iR = 0. 

On the one hand this limiting behaviour can conversely be used in order to label the resonances 
uniquely by two indices m and q, where m G Z and g G N. On the other hand the relations 



(62) are very helpful for finding the resonances numerically. In general, it is quite difficult 
to systematically find all solutions of the equation fm{k) = for finite n in some region 
of the complex plane. For this reason, we start with the resonances for some high value 



of the refractive index n for which the values in (62) serve as good approximations for the 
wavenumbers. We found it convenient to start with n = 50. We then gradually decrease n and 
follow the resonances in the complex plane. We note that the standard routines for computing 
the Bessel and Hankel functions with complex arguments (e.g. NAG, Maple) were not sufficient 
for computing many of those resonances which lie very close to the real axis. We computed the 
results here using the routine of Amos ( |Amo83a] . |Amo83b] ) with the modifications given in 
|YM97] and |YM03] . This is a very sensitive computation for high modal numbers m, and in 
some instances will have false convergences to numbers with very small but positive imaginary 
parts. This happens when the actual solutions are very close to the real axis, and we set the 
imaginary parts to be numerically zero in these cases. Further, we found that the starting value 
'^max should be decreased with increasing m so that by m = 50 we took nmax ~ 7. 

By the symmetry of the Bessel functions we have that km,q = k-rn,q and hence all resonance 
wavenumbers with m 7^ are double degenerate. We denote the multiset of the wavenumbers 
of the inner resonances again by a. Then 

a = {k^^, G C : m G Z, g G N, fm{km,q) = 0} . (63) 



In Fig. 11 we show the resonances for n = 3. One can see that, except for a few of the first 



resonances, they all lie in an infinite strip bounded by the line |DBD"'"08] 

- Imfci? < 7(n), (64) 

where ^ ^ 

li^) = TT 7 , (65) 

2n n — 1 



This bound arises from the asymptotic forms of the Bessel and Hankel functions in (59) for 
large k and m = 0. 

The length spectrum is again given by the Fourier transform F{1) of the density of states 
multiplied with the Gaussian cut-off function W{k) = exp(— tn^fc^). We determined the res- 
onances in a window < Refc < /cmax with fcmax = lOO/ra and chose the cut-off parameter 
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Figure 11: Resonances in the complex plane of the unperturbed circular microcavity with 
refractive index n = 3. Each resonance string corresponds to a fixed radial modal number 
q and the modal number m increases as one moves along the string from the threshold 7(n) 
(dashed line) towards the real axis. 
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t = 10/(nfcinax)^- We hence define 

POO 

F{1) = / d{k) W{k)e-'''^^dk ^ ^e-'"'^"''''-"''^" 
•Jo 



(66) 



m,q 



In contrast to the closed case, the density of states is now given by a sum over Lorentzians, see 



Eq. (11). The integral can be expressed in terms of the error function with complex argument, 



and the approximation in equation (66) holds in our parameter range. 



5.2 Short-wave approximation 

The short-wave approximation for the length spectrum is obtained by inserting the trace formula 



for the density of states into the definition (66). In section 2 it was discussed how the trace 



formula is modified for the open cavity in comparison to the closed cavity. In the smooth part 
of the density of states only the perimeter term changes and we obtain 



ReFo(/) = + 



r(n)L iAl , f il 
+ —erf 
n t \2Vt 



(67) 



where A = nR , L = 2ttR and r(n) is defined in (13) 



The oscillatory part is changed according to equation (15): One has an additional overall 



factor of n and the wavenumber k is replaced by nk. Furthermore, the factors that are due to 
the boundary condition at the cavity wall change. Instead of a Dirichlet phase factor (—1) for 
every reflection we now have to include a Fresnel coefficient for every reflection. For TM modes 
it is given by 
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(68) 



where 9 is the angle of incidence and the critical angle 9c is deflned in (16). This coefficient 



allows for leakage from the cavity and/or entrapment within it (total internal reflection). For 
trajectories in the circular cavity the angle of incidence 9 is the same for each of its reflections. 



and hence we flnd that the oscillatory part for the closed cavity in (23) changes into 
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(69) 



for the open cavity. The index p = (r, w) again labels the periodic orbit families with w 
rotations around the billiard centre and r reflections at the boundary. The quantities in this 



formula have been discussed after the trace formula for the closed cavity (23), except for the 



Fresnel reflection coefficient Rp which is evaluated at the angle of incidence ( 24 ) of each periodic 
orbit. 



We insert the periodic orbit contributions (69) into the definition of the length spectrum 
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and perform the integral analytically. Again we consider only the real part and obtain 



ReF, 



po 



+e" 



V 







2,1 p 








sin 1 




— r< 




/vrr 




sin 




— r 




vY 





X 



(70) 



(Ip + I) 



2t 

{Ip - I) 



+ cos 



Trr 

Y 



{Ip + /) 



2t 



(Txr 
cos y— - r(f 



'2t 
{Ip - I) 



2t 



As mentioned before, the formula becomes less accurate when the angles of incidence of the 
contributing periodic orbits are close to the critical angle |HS02] . We can therefore estimate 
where we expect an error to be located in the length spectrum. If we insert the critical angle 



into the formula for the periodic orbit lengths (25) and use the interpolating formula for the 
number of reflections r = (vr — 26)/2tt of orbits with w = 1 rotations we obtain 
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(71) 



We have that lp{n 



—7- 27iR as n — )■ 1+ and lp{n) — )■ AR as n — > oo, which correspond to the 
perimeter and twice the diameter of the circle, respectively. We therefore have that for suitably 
large refractive index the error is not located near the accumulation point / = 2ttR of the 
periodic orbits but near / = AR at which length we have only the diameter orbit P0(2, 1). 



In Fig. 12 we plot the oscillating part of the semiclassical length spectrum ReFpo(/), together 
with the exact length spectrum minus the smooth part Re[F(/) — Fo(Z)] for three different values 
of the refractive index. Although the periodic orbit structures are well reproduced by the short- 
wave approximation the error is clearly bigger than for the closed cavity, in particular for the 
case n = 1.5. In order to see the error more clearly, we plot in Fig. 13 the difference between 
the exact and semiclassical length spectra. The lengths which correspond to the critical angle 
are I = 5.57 for n = 1.5, / = 4.81 for n = 3, and / = 4.55 for n = 4.5, and we can indeed 
observe large errors at these positions. There are observable errors at other lengths, but they 
correspond to smaller relative errors, because they occur at peaks in the length spectrum in 



Fig. 12 They might be related to the Goos-Hanchen shift [TS02 j whereby the observed length 
at which one has peaks in the exact spectrum is slightly shifted from the length at which one 
has the periodic orbits. We also see an error near / = in the length spectrum which decreases 
with increasing n. This results from only having two terms in the smooth part (12). 



Our test of the trace formula for the dielectric circle goes beyond the results in jBDSOS] 
in that we checked not only the positions of the peaks in the length spectrum but also their 



shapes. Figs. 12 and 13 provide encouraging support for the use of short-wave approximations 
in open dielectric systems and show that the trace formula works as well as can reasonably be 
expected. 



6 The dielectric cavity with point scatterer 
6.1 Exact solution 

As with the closed system we perturb the dielectric cavity by placing a point scatterer at a 
position d = (d, 0) in Cartesian coordinates, see Fig. fl] The exact resonance spectrum for this 
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Figure 12: Length spectra for the circular dielectric cavity, for three values of the refractive 
index n. The exact spectrum is computed from (66) with unperturbed resonances and the 
short-wave approximation from (70) with periodic orbits, (t = 0.001.) 
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Figure 13: Difference between tfie exact lengtfi spectrum of tlie circular dielectric cavity and 
its short-wave approximation for the three values of n in Fig. 12 (t = 0.001.) 
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perturbed open system is again determined by equation (28) where we now have to include the 
regularized Green function of the unperturbed dielectric disk. The regularization of the Green 
function is done exactly as before, by subtracting the logarithmic divergence of the Hankel 
function Hq. The only difference is that its argument contains now nk instead of k. We hence 
obtain that the perturbed resonances, k, are determined as the solutions of f{k) = 0, where / 
is the transcendental function |DMSW08] 

As before, the parameter a is a measure of the strength of the perturbation and replaces the 



and €m in (31). 



coupling constant A appearing in (28). The coefficients and Cm are given in (57) and (58), 



As mentioned in Section |4] one can realize the point scatterer experimentally for small values 
of a by a small disk of radius a with Dirichlet boundary conditions. As discussed in |DMSW09a] 
one can also use a small hole of radius b, filled with dielectric material of refractive index rib. 
These parameters are related to a by 

log — + 7 ~ ,2/2^ 2 ^ • (^3) 

2 b^k^{n-^ — n^) 



The value of a again ranges from to oo, and at these two limiting values the zeros of (72) 



correspond to the unperturbed resonances of Section |5.1[ For finite values of a, one finds 
again that the rank one perturbation changes only one of each pair of degenerate resonances 
for m 7^ of the unperturbed system. This is again the case because the wavefunctions with 



angular dependence sin(m0) in (60) have a nodal line on the x-axis and hence they and their 
corresponding wavenumber are both unaffected by the presence of the scatterer. The other 
resonances have wavefunctions with angular dependence cos(m0) in the unperturbed system, 
and as a changes from to oo their wavenumbers move from one unperturbed value km,q to 
another. We call these the perturbed resonances. The full spectrum a of the perturbed system 
is therefore 

5- = {km,g = km,g e a : -m, g G N} U {Ka,q G C : m G Nq, g G N, f{km,q) = 0} • (74) 

For the labelling we adopt again the rule that km,q has the same m and q values as the unper- 
turbed resonance at which it arrives at a = oo. If we would use the unperturbed resonance 
at a = for the labelling we would again miss the first perturbed resonance. This resonance 
moves from k = fco,i G a at a = oo towards the imaginary axis and then down along it as a is 
decreased. 

The computation of the perturbed resonances is a more delicate procedure than the com- 
putation of the unperturbed resonances in the previous section. In order to find the perturbed 
resonances for a non-zero and finite value of oq we first set a to zero and use an unperturbed 
resonance as an initial starting value in a Newton procedure. (For the first resonance we start 
at a = oo.) Note that we cannot use a bisection method for the open system. Because of the 



logarithmic dependence on a in (72), we then vary a exponentially towards the desired value 
a = oq. The function f{k) has a pole at each unperturbed resonance k = km,q G cr because the 
corresponding in the denominator of its defining sum (72) vanishes at each k = km,q- To 



avoid this numerical problem we scale the function f{k) by the particular Am when computing 
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Figure 14: Level dynamics of the perturbed resonances of the dielectric circular cavity with 
point scatterer {n = 3,d = 0.59) for Rek G [0,4]. The solid (red) lines are the paths of 
the perturbed resonances as the strength varies between a = and a = oo. The dots are 
the unperturbed resonances of the cavity without scatterer and the crosses are the perturbed 
resonances for a = 0.1. For a = 0.1, the position of the first perturbed resonance is not visible 
as its location is deep in the complex plane ai k = 9.5670 x 10~^ — 3. 74351. 



perturbed resonances for the set a. In fact, because a = corresponds to an unperturbed 
resonance, we start the numerical procedure at a small value a = 10~^° and then use an ap- 



proximate solution to ( 72 ) obtained from a perturbative approach. For small (or large) values 



of a this gives the value of the perturbed resonance k^^q as |DMSW0"9aj 

r ^ 7 • ^m{km,q) ^mJ m^km,q^dj^ /'7K^ 

Dra{kra,q) 2R\oga 

where 

D^{k) = {n' - l)JUnkR)HUkR) - . (76) 

As we vary a the perturbed resonances follow paths in the complex plane, each of which 
connects together either two different unperturbed resonances or one unperturbed resonance 
to itself. An example is given in Fig. [14] for low values of k. We also show in the plot the 



unperturbed resonances of Section 5A and the perturbed resonances corresponding to a fixed 
perturbation a = 0.1. We see the first perturbed resonance is moving towards the imaginary 
axis and then much deeper into the complex plane. Some of the paths in Fig. [14] connect 
individual unperturbed resonances to themselves and the others connect different unperturbed 
resonances. There are no instances of a pair of unperturbed resonances being connected in a 
loop-like manner. None of the paths cross (due to level repulsion) but there can be instances 
of two paths becoming very close to each other, and hence a fine discretisation of a is required. 
Furthermore, the choice of the location of the scatterer d determines the overall structure of 
the connections |DMSW09a] . 
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Figure 15: Level dynamics of the perturbed resonances of the dielectric circular cavity with 
point scatterer {n = 3,d = 0.59) for Re/c e [15, 17]. The sohd (red) lines are the paths of the 
perturbed resonances in the set (j\af^ as the strength varies from a = to a = oo and the dashed 
(blue) lines are the level dynamics of the resonances in the set a' which travel to unperturbed 
resonances having small imaginary part as the strength varies from a = oo to a = (see text). 
The dots are the unperturbed resonances of the cavity without scatterer and the crosses are 
the perturbed resonances for a = 0.1. 



Because many of the unperturbed resonances had their imaginary parts set to zero in the 
previous section, the computation of the corresponding perturbed resonances is complicated as 
the quantisation condition (72) is technically unsatisfied at the start of the path. Moreover, 
paths that start at one of these unperturbed resonances near the real axis do not necessarily 
stay in the vicinity of it and hence we cannot arbitrarily take the perturbed resonance to be 
equal to the unperturbed resonance. We verify this by introducing the set of unperturbed 
resonances k = k^+iki with —ki < e. We then omit this set of resonances from the computation, 
vary a from to oo, which is numerically taken to be a = 10^°, and identify the presence of 
any unperturbed resonances k G (j\a^ (that lie deeper in the complex plane) to which no level 
path has reached. The singular nature of the condition (72) at the unperturbed resonances 
prohibits a perturbed resonance from numerically reaching those values. We hence determine 
the 'equality' of two such resonances by using the approximate solution (75). 

We call this set of unperturbed resonances a' and on each element within it perform the 
perturbation computation backwards from a = oo to a = where we always end up at the real 
axis. The unperturbed resonance in which lies closest (up to using (75)) to each termination 
point of each point on the real axis is then removed from a^- The resulting set then contains 
only unperturbed resonances which, when perturbed, follow tiny loops near the real axis back 
to themselves and is hence the set for which we take the perturbed resonance to be equal to 
the unperturbed resonance, regardless of the value of a. The perturbation computation is then 
repeated for a = — t- ao on the set <7\<7e and for a = oo — ag on the set a', the results of which 
contribute, along with the set a^, to the desired set of all perturbed resonances a in (74). 

In Fig. 15 we show the level dynamics of perturbed resonances for k in the range 15 to 17, 
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Figure 16: Exact length spectra for the unperturbed (full line) and perturbed (dashed line) 
circular dielectric cavities, {d = 0.59, a = 0.1, t = 0.001, /imax = 0.) 



and illustrate the above method by highlighting level curves and resonances as computed from 
the sets a\ae or a' separately. As an illustration, for the calculations in the next section we 
will use resonances in the range Ren A; G [0,100] for n = 3, d = 0.59, and a = 0.1. In our 
calculations to find them we used 1278 unperturbed resonances (of which 32 had m = 0), the 
set contained 432 perturbed resonances, the set cr\o"e contained 750 and the set a' contained 
98, with a total of 1280. The discrepancy of 2 comes from an extra perturbed resonance at 
each end of the wavenumber interval. 

Because the diffraction coefficient entering the semiclassical amplitudes changes with the 
additional factor of n, the optimal value of the strength of the scatter a in a given fc-range also 
changes. We found that taking a = 0.1 gives strong diffractive peaks in the length spectra, 
which is smaller than the value a = 1 taken for the corresponding closed system. In our plots 
we use the same scatterer position d = 0.59 as in the closed billiard. 

The exact length spectrum is of the same form as before and is given by 

POO _ 

F{1)= P--/2j(A;)H/(A;)e-'"'=MA;, (77) 
Jo 



where d is the level density of the perturbed system. In Fig. 16 we plot the perturbed length 



spectrum (77) alongside the unperturbed spectrum (66) for a scatterer position d = 0.59, 
scatterer strength a = 0.1 and /^max = 0. As will be discussed in the next section the difference 
can be attributed to diffractive orbits. In order to visualize the signature of diffractive orbits in 
the length spectrum we will again consider in the following the difference, AF{1) = F{1) —F{1), 
between perturbed and unperturbed length spectra . 

6.2 Short-wave approximation 

Similar to the modification of the periodic orbit terms in the trace formula when we moved 
from the closed to the open system, we must also modify the diffractive orbit terms. According 



to equation (15) this consists of adding an overall factor of n to the diffractive part (47) of the 
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Figure 17: Difference between length spectra for perturbed and unperturbed dielectric cavities. 
The exact result (full line) is compared to the short-wave approximation (dashed line) that is 
obtained from diffractive orbits using uniform approximations, {d = 0.59, a = 0.1, /imax = 4, 
t = 0.005.) 



closed system, changing k to nk, and replacing the Dirichlet phase factors by Fresnel reflection 
coefficients. This results in 



ddo{k) = ^ A^exp |i (^nkLd - "y" - ^^^^ | 



+ c.c. 



(78) 



where the amplitudes are given by 
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(79) 



and the sum runs over all diffractive orbits as presented in Section 4^ The refection coefficient 
is now evaluated at the angle of incidence of each contributing diffractive orbit (which is complex 
valued for ghost orbits in a tangent bifurcation). The diffraction and reflection coefficients are 
both complex valued. 

The difference AF{1) between the length spectra of perturbed and unperturbed systems is 
given only in terms of the diffractive orbits in the short-wave approximation, because the peri- 
odic orbit terms are the same in both systems and cancel. Hence we obtain the approximation 



AF(/) 



fc^— /2(ido(A;)Vr(A;)e 



dk. 



(80) 



which has to be evaluated numerically. The quantity /imax is again the chosen maximum number 



of encounters with the scatterer of the diffractive orbits that we enter into the sum. In (80) we 



neglect the small modification of the smooth part of the density of states due to the scatterer. 



In Fig. 17 we plot the length spectrum difference alongside its semiclassical counterpart (80) 



given above. We see that the peaks in AF(/) can clearly be identified with contributions of the 
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Figure 18: Discrepancy between the exact and semiclassical length spectrum differences in 



Fig. 17 {d = 0.59, a = 0.1, /z^a^ = 4, t = 0.005.) 



diffractive orbits. We found it necessary to choose a higher scahng factor t 
previous sections, because we encountered larger fluctuations. 



50/kl 



than in 



As with the perturbed circular billiard, the diffractive contribution ( 78 ) to the trace formula 



breaks down near bifurcations of the diffractive orbits, and we again applied uniform approx- 
imations to deal with these cases. The bifurcation scenarios are the same as for the circular 
billiard and are given in Tab. [T]and we use the same uniform approximations as for the closed 
system, modifying the amplitudes of the orbits involved in each bifurcation appropriately. The 



semiclassical curve in Fig. 17 includes these uniform approximations. We have again included 
the uniform approximation for the tangent bifurcation of the D0(4, 1) and D0(4, 1)' orbits, 
and for the pitchfork bifurcation involving the orbits D0(3, 1) and P0(3, 1). We see in the 
figure that there is a greater discrepancy between the exact and semiclassical length spectra 
corresponding to these orbits than in the analogous figure (Fig. [s]) for the closed billiard. 

In the previous section we saw that the semiclassical approximation to the length spectrum 
of the unperturbed dielectric cavity is poorer at lengths corresponding to periodic orbit whose 
angle of incidence with the boundary is near the critical angle. We gave a formula (71) that 



specifies this length as a function of the refractive index n. We now perform a similar calculation 
for diffractive orbits. If we insert the formula for the critical angle into the formula for the 
diffractive orbit lengths (41), restrict to w = 1 rotations and use an interpolating formula for 



the number of reflections r in terms of the angle of incidence we obtain 
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n 



2 arcsin(l / nd) + tt 
TT - 2er(n) 
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(81) 



n 



which is defined for n > 1/d. At the limiting value n = 1/d it agrees with Ip 
Idin) — )■ 2(1 + (i) as n — )■ oo. For our numerics, the relevant value is /rf(3) ~ 4.30 for d 



in (|71|), and 
0.59. 



At this length we would expect an inaccuracy in the diffractive contribution to the semiclassical 
length spectrum. 



In Fig. 18 we show the discrepancy between the exact and semiclassical length spectra 
differences of Fig. 17 The error at the length that corresponds to the critical angle does not 
play such a big role in this plot. The discrepancy comes mainly from the fact that the peak 
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heights are not completely reproduced by the semiclassical plot, as can be seen in Fig. 17 As 



discussed before, however, there is also a slight uncertainty about the numerical error of some 
of the resonances near the real line. All in all, the semiclassical approximation captures the 
peak structure of the length spectrum quite well. 



7 Conclusions 

The aim of this article was an investigation of trace formulas for dielectric cavities. Such a 
trace formula was proposed recently for TM modes in quasi two-dimensional geometries, and 
it connects the spectrum of the Feshbach or interior resonances in the short-wave limit to the 
ray dynamics inside the cavity [BDSOSj . 

For the dielectric circle we evaluated the Fourier transform of the resonance spectrum, the 
length spectrum, which has peaks at the periodic rays inside the cavity. With more than 
2000 resonances we could probe the short-wave regime. In contrast to previous studies we 
compared not only the position but also the shape of the peaks. We found that the agreement 
with the short-wave approximation is remarkably good, except for expected inaccuracies for 
orbits with reflection angles near the critical angle for total internal reflection. Otherwise, the 
approximation is of similar quality as for the closed cavity. We also saw good agreement for 
orbits which are not in the regime of total internal reflection which were found difficult to see 
in previous articles |BBD+10[ lBDD+09j . 



We then considered a dielectric circular disk with an additional point scatterer. This re- 
quired a generalization of the trace formula which includes additional contributions from diffrac- 
tive rays that start and end at the scatterer. This system allows for a more probing test of 
the trace formula. On the one hand it is a non-integrable system, but one can nevertheless 
determine a large number of resonances by a Green function method. On the other hand the 
contributions of isolated diffractive rays are at least an order of smaller than those of the 
continuous families of periodic orbits in the circle. For example, in |BDD"'"09j it was found 
difficult to get a good quantitative agreement with the trace formula in systems with isolated 
periodic orbits. 

In order to isolate the peaks at diffractive orbits from those at periodic orbits we considered 
the difference of the length spectra for dielectric circular disks with and without scatterer. The 
short-wave approximation for this difference is solely given by diffractive orbits. Our numerical 
results showed that the short-wave approximation works well also in this case. The position 
and shape of the peaks are well reproduced by the trace formula. Only the agreement in 
the height of the peaks was not as good as in the comparable case of the closed cavity with 
a point scatterer. Some diffractive orbits were close to a bifurcations and we improved the 
approximation by using uniform approximations. We also saw clear peaks at diffractive orbits 
which are not in the regime of total internal reflection. 

In summary, we found that trace formulas for dielectric cavities work well in the short-wave 
regime. Even higher order contributions from diffractive rays can clearly be identifled. One 
possible improvement is to include modiflcations of the Fresnel coefficient for reflections near 
the critical angle |BDD"'"09] . As a possible application of our results one could use the resonance 
spectrum to locate the position of a defect within an optical cavity. 
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